// to calculate the single proton knockout cross sections of C16
// @239 MeV/nucleon using MOMDIS and Shell Model
// Created: May 22, 2022 by Yazhou Sun

const double A = 16., Sn = 4.2503, Sp = 22.5531;
const double csExp = 82.7639, dcsExp = 8.3707;

const int n = 1;
const char *orbit[n] = {"0p3/2_3/2-"};
const double Ex[n] = {0.};
const double c2s[n] = {2.95}; // cpc.xu2022
// execute momd to get the following results //
const double csTh[n] = {};




void rate(double N1, double dN1, double N2, double dN2);
void csC16_p(){
	//////////////////////////////////////////////////////////////////////////////////
	double csThTotal = 0., ExAvrg = 0., dSEff = Sn-Sp;
	double sigma[n]{}, N[n]{};
	map<char, int> ll;
	ll['s'] = 0; ll['p'] = 1; ll['d'] = 2; ll['f'] = 3; ll['g'] = 4;
	for(int i = 0; i < n; i++){
		int n, l, j; char lm;
		sscanf(orbit[i], "%d%c%d%*s", &n, &lm, &j);
		if(ll.count(lm)) l = ll.at(lm);
		else{ cout << "ERROR: lm = " << lm << " not found." << endl; getchar(); }
		N[i] = 2.*n+l;
		sigma[i] = pow(A/(A-1.), N[i])*c2s[i] * csTh[i];
		ExAvrg += sigma[i]*Ex[i];
		csThTotal += sigma[i];
	} // end for over i
	ExAvrg /= csThTotal;
	dSEff += ExAvrg;

	cout << fixed << setprecision(6);
	cout << "theoretical c.s.:" << endl;
	cout << left;
	for(int i = 0; i < n; i++) cout << setw(12) << orbit[i] << sigma[i] << setw(10) << " mb" << endl;
	cout << "Total theoretical c.s.: " << csThTotal << " mb" << endl;
	cout << "\033[32;1mRs: ";
	rate(csExp, dcsExp, csThTotal, 0.);
	cout << "\033[0m";
	cout << "dSEff: " << dSEff << endl;
} // end main function

void rate(double N1, double dN1, double N2, double dN2){
  const double rate = N1 / N2;
  const double rel = sqrt(pow(dN1 / N1, 2) + pow(dN2 / N2, 2));
  const double sigma = rate * rel;
  // cout << "relative error: " << rel << endl;
  cout << rate << "(" << sigma << ", " << rel*100. << "%)" << endl;
} // end function rate
